#!/usr/bin/env python
# coding: utf-8

# In[1]:


import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import math
from scipy.stats import skew
from matplotlib.lines import Line2D


# In[2]:


df = pd.read_csv("D:/research/recordbreak/record45/indus_record45_inc.csv")


# In[3]:


# Assuming 'year' is the name of the column containing the years
df['year'] = pd.to_datetime(df['year'], format='%Y').dt.year
df.set_index('year', inplace=True)


# In[4]:


print(df)


# In[5]:


dr_theo = df[df.index <= 2099][['dr_theo']]
print(dr_theo)


# In[7]:


fr_theo = df[df.index <= 2099][['fr_theo']]
print(fr_theo)


# # Increasing

# In[8]:


dr_obs_inc = df[df.index <= 2022][['dr_obs']]
print(dr_obs_inc)


# In[9]:


fr_obs_inc = df[df.index <= 2022][['fr_obs']]
print(fr_obs_inc)


# In[12]:


dr_inc_boot = [f"dr_inc{i}" for i in range(1, 101)]

# Select the specified columns and store them in a new DataFrame
dr_boot_inc = df.loc[df.index <= 2099, dr_inc_boot]
print(dr_boot_inc)


# In[13]:


fr_inc_boot = [f"fr_inc{i}" for i in range(1, 101)]

# Select the specified columns and store them in a new DataFrame
fr_boot_inc = df.loc[df.index <= 2099, fr_inc_boot]
print(fr_boot_inc)


# In[14]:


plt.figure(figsize=(14, 8))
percentile_5_dr = np.percentile(dr_boot_inc, 5, axis=1)
percentile_95_dr = np.percentile(dr_boot_inc, 95, axis=1)
percentile_50_dr = np.percentile(dr_boot_inc, 50, axis=1)

# Plotting the drought record-breaking events decay rate
plt.plot(dr_obs_inc.index, dr_obs_inc.values, label='Observed', color='blue', alpha=0.8, linewidth=3)

for col in dr_boot_inc.columns:
    plt.plot(dr_boot_inc.index, dr_boot_inc[col], color='grey', alpha=0.06)

# Plot 5th and 95th percentiles
plt.plot(dr_boot_inc.index, percentile_50_dr, color='orange', linewidth=3, alpha=0.9, label='RCP 4.5')

plt.plot(dr_boot_inc.index, percentile_5_dr, label='5-95 CI', color='black', linestyle='dashed', alpha=0.7, linewidth=2)
plt.plot(dr_boot_inc.index, percentile_95_dr, color='black', linestyle='dashed', linewidth=2, alpha=0.7)
plt.plot(dr_boot_inc.index, percentile_50_dr, color='orange', linewidth=3, alpha=0.9)

# Theoretical
plt.plot(dr_theo, label='Theoretical (1/n)', color='black', linewidth=2)

# Adding labels and title
plt.xlabel('Water Year (Apr-Mar)', fontsize=18)
plt.ylabel('Record Rate (Log)', fontsize=18)
plt.title('Indus 6-Models Annual-Min increasing Record-breaking DNN Ensemble', fontsize=20)

# Increase axis ticks size
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)

plt.yscale('log')

# Adding a legend
plt.legend(prop={'size': 14})

# Display the plot
plt.show()


# In[15]:


plt.figure(figsize=(14, 8))
percentile_5_fr = np.percentile(fr_boot_inc, 5, axis=1)
percentile_50_fr = np.percentile(fr_boot_inc, 50, axis=1)
percentile_95_fr = np.percentile(fr_boot_inc, 95, axis=1)

# Plotting the flood record-breaking events decay rate
plt.plot(fr_obs_inc.index, fr_obs_inc.values, label='Observed', color='blue', alpha=0.8, linewidth=3)

for col in fr_boot_inc.columns:
    plt.plot(fr_boot_inc.index, fr_boot_inc[col], color='grey', alpha=0.1)

# Plot 5th and 95th percentiles
plt.plot(fr_boot_inc.index, percentile_50_fr, color='orange', linewidth=3, alpha=0.9, label='RCP 4.5')

plt.plot(fr_boot_inc.index, percentile_5_fr, label='5-95 CI', color='black', linestyle='dashed', alpha=0.7, linewidth=2)
plt.plot(fr_boot_inc.index, percentile_95_fr, color='black', linestyle='dashed', linewidth=2, alpha=0.7)

# Theoretical
plt.plot(fr_theo, label='Theoretical(1/n)', color='black', linewidth=2)

# Adding labels and title
plt.xlabel('Water Year (Apr-Mar)', fontsize=18)
plt.ylabel('Record Rate (Log)', fontsize=18)
plt.title('Indus 6-Models Annual-Max increasing Record-breaking DNN Ensemble', fontsize=20)

# Increase axis ticks size
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)

plt.yscale('log')

# Adding a legend
plt.legend(prop={'size': 14})

# Display the plot
plt.show()


# In[16]:


# Calculate the ratio between corresponding percentiles of two datasets
ratio_5 = np.divide(percentile_5_fr, percentile_95_dr)
ratio_50 = np.divide(percentile_50_fr, percentile_50_dr)
ratio_95 = np.divide(percentile_95_fr, percentile_5_dr)


# In[17]:


ratio_obs_inc = df[df.index <= 2022][['ratio_obs_inc']]
print(ratio_obs_inc)


# # RCP 8.5

# In[19]:


df1 =pd.read_csv("D:/research/recordbreak/record85/indus_record85_inc.csv")


# In[20]:


# Assuming 'year' is the name of the column containing the years
df1['year'] = pd.to_datetime(df1['year'], format='%Y').dt.year
df1.set_index('year', inplace=True)
print(df1)


# In[22]:


dr85_boot_inc = [f"dr_inc{i}" for i in range(1, 101)]

# Select the specified columns and store them in a new DataFrame
dr85_boot_inc = df1.loc[df1.index <= 2099, dr85_boot_inc]
print(dr85_boot_inc)


# In[23]:


fr85_boot_inc = [f"fr_inc{i}" for i in range(1, 101)]

# Select the specified columns and store them in a new DataFrame
fr85_boot_inc = df1.loc[df1.index <= 2099, fr85_boot_inc]
print(fr85_boot_inc)


# In[24]:


plt.figure(figsize=(14, 8))
percentile_5_dr85 = np.percentile(dr85_boot_inc, 5, axis=1)
percentile_95_dr85 = np.percentile(dr85_boot_inc, 95, axis=1)
percentile_50_dr85 = np.percentile(dr85_boot_inc, 50, axis=1)

# Plotting the drought record-breaking events increase rate
plt.plot(dr_obs_inc.index, dr_obs_inc.values, label='Observed', color='blue', alpha=0.8, linewidth=3)


for col in dr85_boot_inc.columns:
    plt.plot(dr85_boot_inc.index, dr85_boot_inc[col], color='grey', alpha=0.03)

# Plot 5th and 95th percentiles
plt.plot(dr85_boot_inc.index, percentile_50_dr85, color='orange', linewidth=3, alpha=0.9, label='RCP 4.5')
plt.plot(dr85_boot_inc.index, percentile_5_dr85, label='5-95 CI', color='black', linestyle='dashed', alpha=0.7, linewidth=2)
plt.plot(dr85_boot_inc.index, percentile_95_dr85, color='black', linestyle='dashed', linewidth=2, alpha=0.7)
plt.plot(dr85_boot_inc.index, percentile_50_dr85, color='orange', linewidth=3, alpha=0.9)

# Theoretical
plt.plot(dr_theo, label='Theoretical (1/n)', color='black', linewidth=2)

# Adding labels and title
plt.xlabel('Water Year(Apr-Mar)', fontsize=18)
plt.ylabel('Record Rate (Log)', fontsize=18)
plt.title('Indus 6-Models Annual-Min increasing Record-breaking DNN Ensemble', fontsize=20)

# Increase axis ticks size
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)

plt.yscale('log')

# Adding a legend
plt.legend(prop={'size': 14})

# Display the plot
plt.show()


# In[25]:


plt.figure(figsize=(14, 8))
percentile_5_fr85 = np.percentile(fr85_boot_inc, 5, axis=1)
percentile_95_fr85 = np.percentile(fr85_boot_inc, 95, axis=1)
percentile_50_fr85 = np.percentile(fr85_boot_inc, 50, axis=1)

# Plotting the flood record-breaking events increase rate
plt.plot(fr_obs_inc.index, fr_obs_inc.values, label='Observed', color='blue', alpha=0.8, linewidth=3)


for col in fr85_boot_inc.columns:
    plt.plot(fr85_boot_inc.index, fr85_boot_inc[col], color='grey', alpha=0.03)

# Plot 5th and 95th percentiles
plt.plot(fr85_boot_inc.index, percentile_50_fr85, color='orange', linewidth=3, alpha=0.9, label='RCP 4.5')
plt.plot(fr85_boot_inc.index, percentile_5_fr85, label='5-95 CI', color='black', linestyle='dashed', alpha=0.7, linewidth=2)
plt.plot(fr85_boot_inc.index, percentile_95_fr85, color='black', linestyle='dashed', linewidth=2, alpha=0.7)
plt.plot(fr85_boot_inc.index, percentile_50_fr85, color='orange', linewidth=3, alpha=0.9)

# Theoretical
plt.plot(fr_theo, label='Theoretical (1/n)', color='black', linewidth=2)

# Adding labels and title
plt.xlabel('Water Year(Apr-Mar)', fontsize=18)
plt.ylabel('Record Rate (Log)', fontsize=18)
plt.title('Indus 6-Models Annual-Max increasing Record-breaking DNN Ensemble', fontsize=20)

# Increase axis ticks size
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)

plt.yscale('log')

# Adding a legend
plt.legend(prop={'size': 14})

# Display the plot
plt.show()


# # Ratio

# In[26]:


ratio85_5 = np.divide(percentile_5_fr85, percentile_95_dr85)
ratio85_50 = np.divide(percentile_50_fr85, percentile_50_dr85)
ratio85_95 = np.divide(percentile_95_fr85, percentile_5_dr85)


# In[ ]:





# # Combined

# In[28]:


plt.figure(figsize=(16, 9))
# Plotting the drought record-breaking events increase rate
plt.plot(dr_obs_inc.index, dr_obs_inc.values, label='Observed', color='black', alpha=0.8, linewidth=4)

# 100 predictions from bootstrapping rcp 4.5
#for col in dr_boot_inc.columns:
    #plt.plot(dr_boot_inc.index, dr_boot_inc[col], color='teal', alpha=0.05)

# Plot 5th and 95th percentiles
plt.plot(dr_boot_inc.index, percentile_50_dr, color='teal', linewidth=3, alpha=0.9, label='RCP 4.5')
#plt.plot(dr_boot_inc.index, percentile_5_dr, color='teal', linestyle='dashed', alpha=0.6, linewidth=1.5)
#plt.plot(dr_boot_inc.index, percentile_95_dr, color='teal', linestyle='dashed', alpha=0.6, linewidth=1.5, label='5-95 CI')
plt.fill_between(dr_boot_inc.index, percentile_5_dr, percentile_95_dr, color='teal', alpha=0.2)

# 100 predictions from bootstrapping rcp 8.5
#for col in dr85_boot_inc.columns:
    #plt.plot(dr85_boot_inc.index, dr85_boot_inc[col], color='red', alpha=0.05)

# Plot 5th and 95th percentiles
plt.plot(dr85_boot_inc.index, percentile_50_dr85, color='red', linewidth=3, alpha=0.8, label='RCP8.5')
#plt.plot(dr85_boot_inc.index, percentile_5_dr85, color='red', linestyle='dashed', alpha=0.6, linewidth=1.5)
#plt.plot(dr85_boot_inc.index, percentile_95_dr85, color='red', linestyle='dashed', alpha=0.6, linewidth=1.5, label='5-95 CI')
plt.fill_between(dr_boot_inc.index, percentile_5_dr85, percentile_95_dr85, color='red', alpha=0.15)

#Theoratical
plt.plot(dr_theo, label='Theoretical (1/n)', color='blue', linewidth=3)
# Adding labels and title
plt.xlabel('Water Year(Apr-Mar)', fontsize=22,fontweight='bold') #fontweight='bold'
plt.ylabel('Probabaility (Log)', fontsize=22,fontweight='bold')
#plt.title('Indus wetting Annual-Min Record-breaking', fontsize=20, fontweight='bold')

# Increase axis ticks size
plt.xticks(fontsize=20,fontweight='bold')
plt.yticks(fontsize=18,fontweight='bold')
# Increase the size of x and y tick lines
plt.gca().tick_params(axis='x', which='major', length=8, width=2)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=8, width=2)

plt.yscale('log')
# plt.xticks(dr85_obs.index[::4])
# Set custom y-tick labels for logarithmic scale
yticks_log = [1, 0.5,0.2,0.1,0.05,0.02,0.01]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)
# Adding a legend
plt.legend( prop={'size': 14, 'weight': 'bold'})

numbers_descending = list(range(138, 0, -1))  # Generating a reversed sequence from 138 to 1
ax = plt.gca().twinx()
ax.plot(dr_theo.index, numbers_descending, label='Theoretical (n)', color='green', alpha=0.001, linewidth=2)
ax.invert_yaxis() 
ax.set_ylabel('Return Period (Log)', fontsize=22,fontweight='bold')
ax.yaxis.get_label().set_rotation(270)
ax.yaxis.labelpad = 20

ax.tick_params(axis='y', labelsize=18)
for label in ax.yaxis.get_ticklabels():
    label.set_fontweight('bold')
    

#boundry thickness
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(2)
plt.yscale('log')

yticks_log_secondary = [1, 2,5,10, 20, 50, 100, 150]  # Adjust these values based on your desired logarithmic scale
yticklabels_secondary = [str(val) for val in yticks_log_secondary]  # Convert values to strings for labeling
ax.set_yticks(yticks_log_secondary)
ax.set_yticklabels(yticklabels_secondary)
ax = plt.gca()
ax.tick_params(axis='x', which='major', length=8, width=2)  # Set the desired length, e.g., 10
ax.tick_params(axis='y', which='major', length=8, width=2)  # Set the desired length, e.g., 10

# Display the plot
plt.show()


# In[30]:


#plt.title('Indus wetting Annual-Min Record-breaking', fontsize=20, fontweight='bold')
plt.figure(figsize=(20, 10))
# Plotting the drought record-breaking events decay rate
plt.plot(dr_obs_inc.index, dr_obs_inc.values, label='Observed', color='black', alpha=0.9, linewidth=5)

# Plot 5th and 95th percentiles
plt.plot(dr_boot_inc.index, percentile_50_dr, color='teal', linewidth=4, alpha=0.9, label='RCP 4.5')
plt.fill_between(dr_boot_inc.index, percentile_5_dr, percentile_95_dr, color='teal', alpha=0.2)

plt.plot(dr85_boot_inc.index, percentile_50_dr85, color='red', linewidth=4, alpha=0.8, label='RCP8.5')
plt.fill_between(dr_boot_inc.index, percentile_5_dr85, percentile_95_dr85, color='red', alpha=0.15)

#Theoratical
plt.plot(dr_theo, label='Theo (1/n)', color='blue', linewidth=3)
# Adding labels and title
#plt.xlabel('Water Year(Apr-Mar)', fontsize=22) #fontweight='bold'
#plt.ylabel('Probabaility (Log)', fontsize=22)


# Increase axis ticks size
plt.xticks(fontsize=30) #,fontweight='bold'
plt.yticks(fontsize=30)
# Increase the size of x and y tick lines
plt.gca().tick_params(axis='x', which='major', length=16, width=4)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=16, width=4)

plt.yscale('log')
# plt.xticks(dr85_obs.index[::4])
# Set custom y-tick labels for logarithmic scale
yticks_log = [1, 0.5,0.2,0.1,0.05,0.02,0.01]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)
# Adding a legend
# Adding a legend
#legend = plt.legend(loc='upper right', prop={'size': 30}, bbox_to_anchor=(0.99, 0.98), borderaxespad=0.)
#legend.get_frame().set_linewidth(4)
#legend.get_frame().set_edgecolor('black')

numbers_descending = list(range(138, 0, -1))  # Generating a reversed sequence from 138 to 1
ax = plt.gca().twinx()
ax.plot(dr_theo.index, numbers_descending, label='Theoretical (n)', color='green', alpha=0.001, linewidth=2)
ax.invert_yaxis() 
#ax.set_ylabel('Return Period (Log)', fontsize=22,fontweight='bold')
ax.yaxis.get_label().set_rotation(270)
ax.yaxis.labelpad = 20

ax.tick_params(axis='y', labelsize=30)
#for label in ax.yaxis.get_ticklabels():
    #label.set_fontweight('bold')
    

#boundry thickness
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)
plt.yscale('log')

yticks_log_secondary = [1, 2,5,10, 20, 50, 100, 150]  # Adjust these values based on your desired logarithmic scale
yticklabels_secondary = [str(val) for val in yticks_log_secondary]  # Convert values to strings for labeling
ax.set_yticks(yticks_log_secondary)
ax.set_yticklabels(yticklabels_secondary)
ax = plt.gca()
ax.tick_params(axis='x', which='major', length=16, width=4)  # Set the desired length, e.g., 10
ax.tick_params(axis='y', which='major', length=16, width=4)  # Set the desired length, e.g., 10
plt.xlim(1960, 2100)
# Display the plot
plt.show()


# In[31]:


plt.figure(figsize=(20, 10))
# Plotting the drought record-breaking events decay rate
plt.plot(dr_obs_inc.index, dr_obs_inc.values, label='Observed', color='black', alpha=0.8, linewidth=5)

# Plot 5th and 95th percentiles
plt.plot(dr_boot_inc.index, percentile_50_dr, color='teal', linewidth=4, alpha=0.9, label='RCP 4.5')
plt.fill_between(dr_boot_inc.index, percentile_5_dr, percentile_95_dr, color='teal', alpha=0.2)

plt.plot(dr85_boot_inc.index, percentile_50_dr85, color='red', linewidth=4, alpha=0.8, label='RCP8.5')
plt.fill_between(dr_boot_inc.index, percentile_5_dr85, percentile_95_dr85, color='red', alpha=0.15)

# Theoretical
plt.plot(dr_theo, label='Theo (1/n)', color='blue', linewidth=3)

# Adding labels and title
#plt.xlabel('Water Year (Apr-Mar)', fontsize=22)
#plt.ylabel('Probability (Log)', fontsize=22)

# Increase axis ticks size
plt.xticks(fontsize=30)
plt.yticks(fontsize=30)

# Increase the size of x and y tick lines
plt.gca().tick_params(axis='x', which='major', length=16, width=4)
plt.gca().tick_params(axis='y', which='major', length=16, width=4)

plt.yscale('log')

# Set custom y-tick labels for logarithmic scale
yticks_log = [1, 0.5, 0.2, 0.1, 0.05, 0.02, 0.01]
yticklabels = [str(val) for val in yticks_log]
plt.yticks(yticks_log, yticklabels)

# Adding a legend
#legend = plt.legend(loc='upper right', prop={'size': 30}, bbox_to_anchor=(0.99, 0.98), borderaxespad=0.)
#legend.get_frame().set_linewidth(4)
#legend.get_frame().set_edgecolor('black')

# Set boundary thickness
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)

plt.xlim(1960, 2100)

# Display the plot
plt.show()


# In[32]:


plt.figure(figsize=(16, 9))
# Plotting the flood record-breaking events decay rate
plt.plot(fr_obs_inc.index, fr_obs_inc.values, label='Observed', color='black', alpha=0.8, linewidth=4)

# 100 predictions from bootstrapping rcp 4.5
#for col in fr_boot_inc.columns:
     #plt.plot(fr_boot_inc.index, fr_boot_inc[col], color='teal', alpha=0.03)
# Plot 5th and 95th percentiles
plt.plot(fr_boot_inc.index, percentile_50_fr, color='teal', linewidth=3, alpha=1, label='RCP 4.5')
#plt.plot(fr_boot_inc.index, percentile_5_fr, color='teal', linestyle='dashed', alpha=0.6, linewidth=1.5, label='5-95 CI')
#plt.plot(fr_boot_inc.index, percentile_95_fr, color='teal', linestyle='dashed', alpha=0.6, linewidth=1.5)
plt.fill_between(fr_boot_inc.index, percentile_5_fr, percentile_95_fr, color='teal', alpha=0.2)

# 100 predictions from bootstrapping rcp 8.5
#for col in fr85_boot_inc.columns:
    #plt.plot(fr85_boot_inc.index, fr85_boot_inc[col], color='red', alpha=0.025)

# Plot 5th and 95th percentiles
plt.plot(fr85_boot_inc.index, percentile_50_fr85, color='red', linewidth=3, alpha=0.75, label='RCP8.5')
#plt.plot(fr85_boot_inc.index, percentile_5_fr85, color='red', linestyle='dashed', alpha=0.6, linewidth=1.5,label='5-95 CI')
#plt.plot(fr85_boot_inc.index, percentile_95_fr85, color='red', linestyle='dashed', alpha=0.6, linewidth=1.5)
plt.fill_between(fr_boot_inc.index, percentile_5_fr85, percentile_95_fr85, color='red', alpha=0.15)

#Theoratical
plt.plot(dr_theo, label='Theoretical (1/n)', color='blue', linewidth=3)
# Adding labels and title
plt.xlabel('Water Year(Apr-Mar)', fontsize=22,fontweight='bold') #fontweight='bold'
plt.ylabel('Probabaility (Log)', fontsize=22,fontweight='bold')
#plt.title('Indus wetting Annual-Max Record-breaking', fontsize=20, fontweight='bold')

# Increase axis ticks size
plt.xticks(fontsize=20,fontweight='bold')
plt.yticks(fontsize=18,fontweight='bold')
# Increase the size of x and y tick lines
plt.gca().tick_params(axis='x', which='major', length=8, width=2)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=8, width=2)

plt.yscale('log')
# plt.xticks(dr85_obs.index[::4])
# Set custom y-tick labels for logarithmic scale
yticks_log = [1, 0.5,0.2,0.1,0.05,0.02,0.01]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)
# Adding a legend
plt.legend( prop={'size': 14, 'weight': 'bold'})

numbers_descending = list(range(138, 0, -1))  # Generating a reversed sequence from 138 to 1
ax = plt.gca().twinx()
ax.plot(dr_theo.index, numbers_descending, label='Theoretical (n)', color='green', alpha=0.001, linewidth=2)
ax.invert_yaxis() 
ax.set_ylabel('Return Period (Log)', fontsize=22,fontweight='bold')
ax.yaxis.get_label().set_rotation(270)
ax.yaxis.labelpad = 20

ax.tick_params(axis='y', labelsize=18)
for label in ax.yaxis.get_ticklabels():
    label.set_fontweight('bold')
    

#boundry thickness
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(2)
plt.yscale('log')

yticks_log_secondary = [1, 2,5,10, 20, 50, 100, 150]  # Adjust these values based on your desired logarithmic scale
yticklabels_secondary = [str(val) for val in yticks_log_secondary]  # Convert values to strings for labeling
ax.set_yticks(yticks_log_secondary)
ax.set_yticklabels(yticklabels_secondary)
ax = plt.gca()
ax.tick_params(axis='x', which='major', length=8, width=2)  # Set the desired length, e.g., 10
ax.tick_params(axis='y', which='major', length=8, width=2)  # Set the desired length, e.g., 10

# Display the plot
plt.show()


# In[52]:


#plt.title('Indus drying Annual-Max Record-breaking', fontsize=20, fontweight='bold')
plt.figure(figsize=(20, 10))
# Plotting the drought record-breaking events decay rate
plt.plot(fr_obs_inc.index, fr_obs_inc.values, label='Observed', color='black', alpha=0.9, linewidth=5)

# Plot 5th and 95th percentiles
plt.plot(fr_boot_inc.index, percentile_50_fr, color='teal', linewidth=4, alpha=0.9, label='RCP 4.5')
plt.fill_between(fr_boot_inc.index, percentile_5_fr, percentile_95_fr, color='teal', alpha=0.2)

plt.plot(fr85_boot_inc.index, percentile_50_fr85, color='red', linewidth=4, alpha=0.8, label='RCP8.5')
plt.fill_between(fr_boot_inc.index, percentile_5_fr85, percentile_95_fr85, color='red', alpha=0.15)

#Theoratical
plt.plot(fr_theo, label='Theo (1/n)', color='blue', linewidth=3, alpha=0.8)
# Adding labels and title
#plt.xlabel('Water Year(Apr-Mar)', fontsize=22) #fontweight='bold'
#plt.ylabel('Probabaility (Log)', fontsize=22)


# Increase axis ticks size
plt.xticks(fontsize=30) #,fontweight='bold'
plt.yticks(fontsize=30)
# Increase the size of x and y tick lines
plt.gca().tick_params(axis='x', which='major', length=16, width=4)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=16, width=4)

plt.yscale('log')
# plt.xticks(dr85_obs.index[::4])
# Set custom y-tick labels for logarithmic scale
yticks_log = [1, 0.5,0.2,0.1,0.05,0.02,0.01]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)
# Adding a legend
# Adding a legend
#legend = plt.legend(loc='upper right', prop={'size': 30}, bbox_to_anchor=(0.99, 0.98), borderaxespad=0.)
#legend.get_frame().set_linewidth(4)
#legend.get_frame().set_edgecolor('black')

numbers_descending = list(range(138, 0, -1))  # Generating a reversed sequence from 138 to 1
ax = plt.gca().twinx()
ax.plot(dr_theo.index, numbers_descending, label='Theoretical (n)', color='green', alpha=0.001, linewidth=2)
ax.invert_yaxis() 
#ax.set_ylabel('Return Period (Log)', fontsize=22,fontweight='bold')
ax.yaxis.get_label().set_rotation(270)
ax.yaxis.labelpad = 20

ax.tick_params(axis='y', labelsize=30)
#for label in ax.yaxis.get_ticklabels():
    #label.set_fontweight('bold')
    

#boundry thickness
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)
plt.yscale('log')

yticks_log_secondary = [1, 2,5,10, 20, 50, 100, 150]  # Adjust these values based on your desired logarithmic scale
yticklabels_secondary = [str(val) for val in yticks_log_secondary]  # Convert values to strings for labeling
ax.set_yticks(yticks_log_secondary)
ax.set_yticklabels(yticklabels_secondary)
ax = plt.gca()
ax.tick_params(axis='x', which='major', length=16, width=4)  # Set the desired length, e.g., 10
ax.tick_params(axis='y', which='major', length=16, width=4)  # Set the desired length, e.g., 10
plt.xlim(1960, 2100)
# Display the plot
plt.show()


# In[35]:


plt.figure(figsize=(20, 10))
# Plotting the drought record-breaking events decay rate
plt.plot(fr_obs_inc.index, fr_obs_inc.values, label='Observed', color='black', alpha=0.8, linewidth=5)

# Plot 5th and 95th percentiles
plt.plot(fr_boot_inc.index, percentile_50_fr, color='teal', linewidth=4, alpha=0.9, label='RCP 4.5')
plt.fill_between(fr_boot_inc.index, percentile_5_fr, percentile_95_fr, color='teal', alpha=0.2)

plt.plot(fr85_boot_inc.index, percentile_50_fr85, color='red', linewidth=4, alpha=0.8, label='RCP8.5')
plt.fill_between(fr_boot_inc.index, percentile_5_fr85, percentile_95_fr85, color='red', alpha=0.15)

# Theoretical
plt.plot(fr_theo, label='Theo (1/n)', color='blue', linewidth=3)

# Adding labels and title
#plt.xlabel('Water Year (Apr-Mar)', fontsize=22)
#plt.ylabel('Probability (Log)', fontsize=22)

# Increase axis ticks size
plt.xticks(fontsize=30)
plt.yticks(fontsize=30)

# Increase the size of x and y tick lines
plt.gca().tick_params(axis='x', which='major', length=16, width=4)
plt.gca().tick_params(axis='y', which='major', length=16, width=4)

plt.yscale('log')

# Set custom y-tick labels for logarithmic scale
yticks_log = [1, 0.5, 0.2, 0.1, 0.05, 0.02, 0.01]
yticklabels = [str(val) for val in yticks_log]
plt.yticks(yticks_log, yticklabels)

# Adding a legend
legend = plt.legend(loc='upper right', prop={'size': 40}, bbox_to_anchor=(0.99, 0.98), borderaxespad=0.,ncol=2)
legend.get_frame().set_linewidth(4)
legend.get_frame().set_edgecolor('black')

# Selecting and adjusting the legend entries for RCP 4.5 and RCP 8.5
#handles, labels = ax.get_legend_handles_labels()
#obs_legend = Line2D([0], [0], marker='s', color='w', markerfacecolor='black', markersize=24)
#ax.legend([obs_legend, handles[0], handles[3]], ['Observed', 'RCP 4.5', 'RCP 8.5'], loc='upper left', prop={'size': 36})
#legend = ax.legend([obs_legend, handles[0], handles[3]], ['Observed', 'RCP 4.5', 'RCP 8.5'], loc='upper left', prop={'size': 36},ncol=3)
#legend.get_frame().set_linewidth(4)
#legend.get_frame().set_edgecolor('black')


# Set boundary thickness
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)

plt.xlim(1960, 2100)

# Display the plot
plt.show()


# In[51]:


plt.figure(figsize=(20, 10))
# Plotting the drought record-breaking events decay rate
plt.plot(fr_obs_inc.index, fr_obs_inc.values, label='Observed', color='black', alpha=0.9, linewidth=5)
# Theoretical
plt.plot(fr_theo, label='1/n', color='blue', linewidth=3, alpha=0.8)

# Plot 5th and 95th percentiles
plt.plot(fr_boot_inc.index, percentile_50_fr, color='teal', linewidth=4, alpha=0.9, label='RCP 4.5')
plt.fill_between(fr_boot_inc.index, percentile_5_fr, percentile_95_fr, color='teal', alpha=0.2)

plt.plot(fr85_boot_inc.index, percentile_50_fr85, color='red', linewidth=4, alpha=0.8, label='RCP 8.5')
plt.fill_between(fr_boot_inc.index, percentile_5_fr85, percentile_95_fr85, color='red', alpha=0.15)



# Adding labels and title
#plt.xlabel('Water Year (Apr-Mar)', fontsize=22)
#plt.ylabel('Probability (Log)', fontsize=22)

# Increase axis ticks size
plt.xticks(fontsize=30)
plt.yticks(fontsize=30)

# Increase the size of x and y tick lines
plt.gca().tick_params(axis='x', which='major', length=16, width=4)
plt.gca().tick_params(axis='y', which='major', length=16, width=4)

plt.yscale('log')

# Set custom y-tick labels for logarithmic scale
yticks_log = [1, 0.5, 0.2, 0.1, 0.05, 0.02, 0.01]
yticklabels = [str(val) for val in yticks_log]
plt.yticks(yticks_log, yticklabels)

# Adding a legend with a single line containing three columns
#legend = plt.legend(loc='upper right', prop={'size': 40}, bbox_to_anchor=(0.99, 0.98), borderaxespad=0., ncol=2)
#legend.get_frame().set_linewidth(4)
#legend.get_frame().set_edgecolor('black')

# Set boundary thickness
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)

plt.xlim(1960, 2100)

# Display the plot
plt.show()


# In[50]:


plt.figure(figsize=(20, 10))

# Plotting the flood record-breaking events decay rate
plt.plot(ratio_obs_inc.index, ratio_obs_inc.values, label='Observed', color='black', alpha=0.9, linewidth=5)

#label='Theoretical=1',
plt.axhline(y=1, color='blue', linestyle='solid',  linewidth=3, label='1/n') #label='Theoretical=1',


plt.plot(fr_boot_inc.index, ratio_50, color='teal', linewidth=4, alpha=0.9, label='RCP 4.5')
#plt.plot(fr_boot_inc.index, ratio_5, color='teal', linestyle='dashed', alpha=0.6, linewidth=1.5, label='5-95 CI')
#plt.plot(fr_boot_inc.index, ratio_95, color='teal', linestyle='dashed', alpha=0.6, linewidth=1.5)
plt.fill_between(fr_boot_inc.index, ratio_5, ratio_95, color='teal', alpha=0.2)

plt.plot(fr85_boot_inc.index, ratio85_50, color='red', linewidth=4, alpha=0.8, label='RCP 8.5')
#plt.plot(fr85_boot_inc.index, ratio85_5, color='red', linestyle='dashed', alpha=0.6, linewidth=1.5, label='5-95 CI')
#plt.plot(fr85_boot_inc.index, ratio85_95, color='red', linestyle='dashed', alpha=0.6, linewidth=1.5)
plt.fill_between(fr85_boot_inc.index, ratio85_5, ratio85_95, color='red', alpha=0.15)


#plt.xlabel('Water Year(Apr-Mar)', fontsize=30) #fontweight='bold'
#plt.ylabel('Probabaility Ratio (Log)', fontsize=30)
#plt.title('Indus drying Annual-Min Record-breaking', fontsize=20, fontweight='bold')

# Increase axis ticks size
plt.xticks(fontsize=36) #,fontweight='bold'
plt.yticks(fontsize=36) #,fontweight='bold'
# Increase the size of x and y tick lines
plt.gca().tick_params(axis='x', which='major', length=18, width=4)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=18, width=4)

# Set logarithmic scale for y-axis
plt.yscale('log')
# Set custom y-tick labels for logarithmic scale
yticks_log = [0.2, 0.5, 1, 2, 5, 10]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)

# Increase the thickness of the outer boundary
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)

# Adding a legend
legend = plt.legend(loc='upper right', prop={'size': 34}, bbox_to_anchor=(0.99, 0.98), borderaxespad=0., ncol=4)
legend.get_frame().set_linewidth(4)
legend.get_frame().set_edgecolor('black')

plt.ylim(0.19, 10)
plt.xlim(1960, 2100)
# Display the plot
plt.show()


# In[ ]:





# In[ ]:





# In[ ]:





# In[ ]:





# In[ ]:




